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The explicit form of perturbation equation for the ^4 Weyl scalar, containing the matter source 
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the Schwarzschild spacetime using in-going penetrating coordinates. As a practical application, we 
focused on the emission of gravitational waves when a black hole is perturbed by a surrounding 
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allow, by means of a spherical harmonic decomposition, to study the problem by means of a one 
dimensional numerical code. 
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I. INTRODUCTION 

One of the problems which has received more attention in Numerical Relativity from early nineties, is the one 
related to the sources of gravitational waves, which in turn could be observed by the gravitational wave detectors. 
The strongest candidates for the current detectors are the strongly gravitating binary systems, involving highly non- 
linear evolutions of either black holes or neutron stars. There are other cases where the geometry of the spacetime is 
determined mainly by a big massive compact body which is perturbed by a smaller source. Typical examples are the 
Extreme Mass Ratio Inspirals (EMRIs), with small compact objects (ie, black holes or stars) falling into a supcrmassive 
black hole, or the accretion of matter unto a black hole, where the black hole is surrounded by a much less massive 
cloud of matter. The changes in the spacetime by these less massive objects (ie, either smaller black holes or the 
accreted disk) can be neglected and the reaction of the central object to the motion of the matter are well described 
as perturbations in that background. The perturbativc methods for determining the generation of gravitational waves 
are expected to correctly describe all the stages of the process. Such methods, together with a detail analysis of the 
sources of the perturbations, are then a promising tool for obtaining the features of the gravitational wave, with a 
much less need of computational resources than those needed in the full numerical evolutions. 

Our final goal is to accurately describe the gravitational waves generated by the motion of either small compact 
bodies (EMRI) or disks of matter in the background of a rotating black hole using the perturbation theory. Specifically, 
we work with the curvature perturbations within the null tetrad formulations developed by Newman and Penrose 
[l|. Basically, some scalar quantities are constructed by means of projections of the Weyl tensor onto a basis tetrad 
formed out from null vectors. Projecting the Bianchi identities on the same null tetrad, we derive the field equations 
for those scalar quantities, and then study the perturbed form of these equations in order to obtain an evolution 
equation for the perturbed scalar quantity. Usually this scalar quantity is the perturbed iff 4 scalar which, as it can 
be inferred from the "peeling off theorem" 0, 0|, is the one which describes the outgoing gravitational wave. 

Since the first derivation of the perturbation equation, there have been several specific works in the same line. 
Indeed in [j| , the authors studied some of the specific features of the gravitational radiation produced by a boosted 
dust shell surrounding a fixed non-rotating black hole. For the cases considered they needed to use a 2d code and the 
spacetime was described with the Boyer-Lindquist coordinates, which introduced several numerical difficulties due to 
the singular slicing. Other works, as for instance @, Q, studied the perturbation equation in the case of a rotating 
black hole, using the so called penetrating coordinates (also known as Kerr-Schild or Eddington-Finkelstcin), avoiding 
the divergence problems of the Boyer-Lindquist description, but the authors worked only in the sourceless case. The 
use of these penetrating coordinates within the perturbation equation including sources was started in , where the 
authors explicitly derived the equations for several choices of the null tetrad. 

In the present work, using the ideas mentioned above, we continue with the description of the gravitational waves 
including a treatment of the sources which generate them. We rederived the perturbation equations for the gravita- 
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tional wave, finding a dependence in that equation on the choice of signature. We concentrate on the non rotating 
black hole case, with dust infalling unto the black hole. Decomposing the dust density in spherical harmonics, we are 
able to deal with non spherical distributions within a spherical description. We perform the numerical evolution for 
different initial distributions of the dust and obtain the corresponding gravitational wave, remarking how these waves 
carry the information of the black hole parameters, in this case the mass, as well as information about the distribution 
of the source which generated them. 

The work is organized as follows: in section HH we present a derivation of the perturbation equation for the 
perturbed "J 4 Wcyl scalar, including the source terms, for a general type D vacuum background spacetime (the black 
hole spacctimcs belong to this type). Then, in lllli we use this equation for the static spherical symmetric spacetimes, 
present the Schwarzschild black hole described in Boyer-Lindquist coordinates as well as in penetrating ones, and 
explicitly derived the evolution equation for $ = r$ 4 including sources. In section IIVI we consider that the source 
term is described by radially infalling dust, make an harmonic decomposition of density, and derive the corresponding 
equations. We also decompose $ in terms of spin weighted spherical harmonics obtaining a clear correspondence 
between the modes of matter and the gravitational ones. The perturbation equation can be separated in radial- 
temporal and angular parts, resulting in a radial-temporal equation for each mode, for which we construct a first 
order system of equations. In section [V] we describe the numerical code used to evolve this system of equations, we 
first present a well known case to show that the code is working properly, we also analyze the infall of one gaussian 
pulse and the gravitational response due to the variation in the width of such pulse and we also present the case 
of three infall pulses of matter and the case of two consecutive pulses varying the initial separation between them. 
Finally in section IVT1 we present a discussion of the results obtained. 



II. PERTURBATION EQUATION 

The spinor formulation was introduced by Newman and Penrose [l[ and used by Teukolsky in [J] to derive a master 
equation for different fields in a Kerr background. One of the central ideas of such formulation consists in choosing 
a tetrad of null vectors Z a ^. These null vectors allow to define directional operators (ie, as the projections of the 
partial derivative along each null vector) and the covariant derivatives of such null vectors projected on themselves, 
which plays the role of the Christoffel symbols within this formulation. For a complete review on the subject and its 
properties, the reader is referred to [9I. [To|. 

Usually the tetrad is defined with the following choice for the the four null vectors: two along the light cone and 
the other two in the perpendicular plane to the cone. These last two null vectors are usually defined in terms of a 
complex one and its complex conjugate. The two null vectors along the light cone are taken as real quantities and 
defined such that they are oriented towards the future and, in the asymptotic region of an hypersurface of constant 
time, one points outward and the other inward the hypersurface. As it can be seen, there is a large room for several 
definitions, which might raise some confusion. We make the following choice: we label the two real vectors along the 
light cone as — I 1 *, Z]/ 1 = and choose them both future directed and with pointing outward, and fc M pointing 
inward in the asymptotic region of an hypersurface of constant time. This choice of orientation is consistent with the 
results obtained for the behavior of the Weyl scalar in the asymptotic regions mentioned in the introduction. The 
other two null vectors, in the perpendicular plane, will be denoted by a complex vector and its conjugate m*^, so 
that Z<f = m» and Z^ = m***. 

After these choices, there is still the issue of the normalization of the products l^ 1 k^, and m M m* M . The rest of the 
products are zero by construction. All the properties mentioned above can be expressed in the following equations: 

Za!* Z bfi = 7] ab , g^y = rf b Z a ^Z bv ), (1) 

with rjab a matrix with the form: 



1]ab 



/ 77 

r] 

-77 

V -t] 



(2) 



with 77 a constant related to the normalization. It is usual to set this constant to one, as it was done in fToj ] , and 
taken from there in Q to derive the perturbation equation that we mentioned. However, this choice depends on the 
signature of the spacetime. Indeed, if the spacetime has signature (+,—,—,—), as was the case in the works just 
mentioned, then this choice is consistcn with Eqs. |T]). However, if the spacetime is described using the signature 
(—,+,+,+), as it is common nowadays, we have to choose rj = —1, in order to remain consistent with Eqs. JT]). Some 



discussion on the change of the tetrad due to the signature can be found in 11 1 
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We are using the definitions for the directional operators, spinor coefficients, and projections of the Wcyl tensor as 
given in lief . In this way, for the directional operators we have: 

D = l"0„, A = k^ dfj_, £ = m"d M , (3) 

while the spinor coefficients can be written as: 

K s = m^l Ku r, t s = m fl l/i-vk", a s = m* 1 Z M; „ to" , p s = m M to*", 

tt s = fc M m* WI/ r, v a — k^ to* KV k v ', fj, s = k 1 * m* w „ m u , A s = k^ m* M; „ m* v , 

e s = \ (fc M l mv + ro" m\., u ) l", ls = i (*" i w „ + to" mV) 

& = - + m tl m\. v ) to", a s = - (k» l KV + to' 1 m\, v ) to*", (4) 

where we have added a subindcx s to avoid confusion with other symbols which will be defined later in the work. For 
the projections of the Weyl tensor C^xt, the following five complex quantities, known as Weyl scalars are defined: 



*0 = 




Fm l/ l x m T , 


*1 = 




l"k v l x m T , 


*2 = 




rmm k , 


*3 = 




F k v m* X k T , 


*4 = 




k^m* 1 - 'k x m* 



(5) 

In order to obtain the perturbed equations, one starts from the Bianchi identities and the definition of the Ricmann 
tensor, namely 

These equations can be projected on the null tetrad, leaving unspecified the normalization constant r]. With the 
above definition (|4|5j) . the equations ((6|) translate in equations for the Weyl scalars and spinor coefficients Q. The 
perturbed expression of these equations is computed for the case of vacuum type D spacetimes to obtain a master 
equation for the perturbed Weyl scalar ^4 including the source terms Q: 

[(A + 7] (4/i a + + 3 7s - 7.*)) ( D - V (Ps -4e.)) - {&* + V (3a s + P.* +4tt s - r s *)) (5 + v (4/3, - r s )) 
-3 7?* 2 ]*4 (1) =V§T 4 . (7) 

with K = 8 7r in geometrized units. The source terms, denoted by T4, are given as follows: 

4 — 1 J-kk + l J-km* + 1 J-m*m*i \°l 

where we have introduced the following operators T ab acting on the projections of the perturbed source term v on 
the null tetrad, 

f kk = -(5* +7] (3« s +ft'+4^-r/)) (5*+ V (2a s + 2f3 s *-T s *)), 
f km * = (A + r? ( 4/is + Ms * + 3 7s _ 7s * )) ( 6 *+2r) (a s -r s *)) + 
(5* + V (3a s +P s * +4tt s -t s *)) (A + 2rj (/i s *+ 7s )), 
f m ' m ' = -(A + jj (4 Ms+/ x s *+3 7s -7 s *)) (A + r, (ii s * +2 ls -2 ls *)). (9) 

By choosing 77 = 1, we recover the perturbation equation for ^4 derived by Teukolsky [J by using the signature 
(+, — , — , — ). However, for the other choice of signature, we must take r/ = — f in order to have a consistent description 
of the gravitational perturbation. We remark that the changes in the perturbation equation due to the choice of 
signature are very noticeable and important at this level of the perturbation equation. However, for a specific case (ie, 
for a given spacetime with a given selection of the signature) , the final perturbation equation in a explicit coordinate 
system is invariant with respect to the choice of signature. It is only at the level of the perturbation equation within 
the null tetrad formulation, that care must be taken with respect of the election of signature made. 
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III. SPHERICAL SYMMETRIC SPACETIMES 



In this section we study the perturbation equation ([7]), with the sources given by eqs. (|8l9p . for a static spherically 
symmetric space time. In this case, the most general line clement can be written as: 

ds 2 = - (a 2 - 7 2 f3 2 ) dt 2 + 2 7 2 /3 dt dr + 7 2 dr 2 + r 2 dtt 2 , (10) 

where dfl 2 = d6 2 + sin 2 9 dip 2 is the solid angle element, and the lapse, a, the radial component of the shift vector, /?, 
and the rr-component of the metric, 7 2 , are functions of the r coordinate only. Notice that we have already chosen 
our signature, so that in the perturbation equation we have to use rj = — 1. 

The next step is to construct the null tetrad for this spacetime by using eq. ([1]) with r\ = — 1. With respect to the 
two real null vectors, we choose them to describe the temporal-radial part of the spacetime, so that they only have 
t and r components. Wc define both of them to point to the future in the asymptotic flat region, and in a timclikc 
hypersurface Z** will point outward while A: M will point inwards. Under these conditions, there is left only one free 
function, which we denote by ko = ko(r). The remaining two null vectors have now to describe the angular part. After 
the normalization conditions are imposed on them, there is left another free function which describes a rotation in the 
angular plane. Wc have chooscn this function to be equal to one, since other choices have no simplifying consequences 
on the final perturbation equation. In this way, our general null tetrad for the line element (|10[) reads as: 

'^T-JiTl 1 '--/ 3 ' ' ) • ^ = *b(l,---A0,0) , m" = -^ r (O,O,l,icsc0). (11) 
Within this tetrad, the only non-zero spinor coefficients and Weyl scalars, given by eq. ([3]), are 



Ps = o Tl~ > V 



fco(f+/?) 



d r (ra^fp s ) 



2 r or fco r 2 a 7 

_ d r (ra ll i s ) _ cote _ _ d r (r 2 a 1Ps ii s ) 

Is ^ i Ps — - 7= , a> s — ~p s , W 2 — . (12) 

2a-f 2\/2r raj 

As mentioned in the introduction, within the spherical symmetric spacetimes, the perturbation equation |(7J) has been 
studied to describe the gravitational waves of a non rotating black hole. Among these studies, there are those which 
have used the Boycr Lindquist coordinate system including sources, and those using penetrating coordinates in the 
sourceless case. Let us give a brief review of these works. 

In Q , the authors made pioneering numerical studies of the matter perturbing a Schwarzschild black hole described 
in the Boyer-Lindquist coordinates, which are obtained from the general line element (1101) . with 



M 
r 



a 2 = 1-2 — , 7 2 = — , = ■ (13) 



In that work they used the Kinnersley tetrad, which can be obtained from our general tetrad form (|11[) with the 
choice of metric coefficients (|13|) and ko = \- A direct substitution of these metric coefficients and tetrad vectors on 
Eqs. (fT2l [7)). allows us to obtain the spinor coefficients and then the following perturbation equation: 

+ r(r-2M)— + 4 — i '- — + 2 (3 r - 7 M) — +4 1-4 



r - 2 M dt 2 v ' dr 2 r - 2 M dt v ' dr 

d 2 Id 2 n d ,.cos0 d cos 2 6 + 1 



cot 6 — - Ai- 



de 2 sin 2 6 dp 2 d6 sin 2 9 d V «- 2 



sm 



* 4 (1) = 2Kr 2 T 4 . (14) 



This equation is equivalent in the non-rotating case to the one used in Q , as well as the one obtained by Teukolsky [i[ , 
for the non-rotating case. We do not follow such derivations because to our purpose this equation is enough. Indeed, 
we confirm that there is a divergence problem at the horizon (ie, r = 2M), as expected from the Boyer-Lindquist 
coordinates. This coordinate singularity forces the introduction of the tortoise coordinate in order to overcome this 
divergence. Also, there appears a remarkable problem outside the horizon at r = 3M; there is a change of sign in 
the coefficient of the first temporal derivative. The authors in Q needed to include dissipation techniques to prevent 
the appearance of numerical instabilities at this point. The coefficients of the source terms also present these kind 
of coordinates problems and was needed the ingenuity of the authors to perform the numerical evolutions and to be 
able to study the gravitational answer of the black hole to the falling of several cases of dust shells studied in that 
work. Wc consider interesting to remark that the signature chosen in these two works 0, H| was (+,—,—,—) . As 
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mentioned above, when the derivation is done consistently, the final perturbation equation (ie, the cq. (|14[) in this 
case), is independent of the choice of signature. 

There are also works where the perturbation equation has been numerically solved even for the case of a rotating 
black hole, either by using Boyer-Lindquist coordinates @ and facing similar problems as those mentioned above, or 
by using a coordinate system where the metric coefficients are free of coordinate singularities in the radial sector 0] • 
In this case, the perturbation equation is also free of coordinates singularities and the inner boundary can be set to be 
inside the black hole horizon, allowing the use of the excision technique to avoid dealing with the physical singularity. 
This description results in a system of equations whose numerical implementation is much simpler that when using 
Boyer-Lindquist coordinates. Both of these works considered the vacuum spacetimes, studying the evolution of given 
initial configurations of the perturbed Weyl scalar 

It seems clear that the best option is to continue these works by considering realistic sources while using the coordi- 
nate description free of singularities . In the rest of the present work we perform such program for the Schwarzschild 
case, working along the lines initiated by Papadopoulos and Font p|, but in penetrating coordinates. The line element 
of the Schwarzschild black hole, written by using the ingoing Kerr-Schild coordinates, has the form: 

, , / 2M\ , 9 4M , , / 2M\ , 99 ,^, 

ds 2 = - 1 dt 2 + dt dr + 1 + dr 2 +r 2 dQ 2 . (15) 

\ r J r \ r J 

where, by comparison with the general line element (jlOp . we see that the lapse, the radial component of the shift 
function and the rr-componcnt of the metric are given by 



2M n 2M , 1 

7 2 = 1 + , = — , a 2 = — , 16 

t v 'y 



so all of them are regular outside the singularity. By inspecting the expression (fTTj) for the tetrad in the general 
spheric case becomes clear that fc = 1 gives a particular simple description. This choice differs with the one chosen in 
J7j , which for the Schwarzschild case reduces to consider fc = . The explicit form of our tetrad in these coordinates 
is given by: 

1 / 2M 2AI \ 1 
i" = - M +—,l-—0 J 0j , fc^ = (1,-1,0,0) , m" = — (0,0,1,* esc 0) , (17) 

as well as the non-zero spinor coefficients and Weyl scalar: 

1 r - 2 M M 



H a = - , p s = 



r 2r 2r 

cot 6 M , . 

2 V 2 r r d 

A straightforward substitution of these quantities in the perturbation equation leads to the final equation for the 
perturbed scalar of ^4: 

[□* +n 64> ] ^ 4 (1) ^2Kr 2 T i . (19) 

where 

□* = ~(r 2 + 2Mr) — + ( r 2 - 2 M r) _ + 4Mr — + 2 (2 r + 3 M) - + 6 (r - M) - + 4 (20) 

d 2 1 d 2 d cos9 d l + cos 2 6i , N 

□e v = ^ + -^T^ + cot0 M -4i— — Dir -2 — rT — . 21 
06 2 sin 9 dip 2 89 sm z 6 dip snr 6 

On the other side, the explicit form of the operators defined by the eqs. (|8I9[) for the source term T4 take the form: 

t kk = -3^2 3-1 §o , (22) 
*"--4(±-i-l)*-l. P3) 



dt dr 

T """ - -l'^-»ra;+^-;i=- = l + 3l . < 2 ") 



(-- 


~) 




\dt 


dr J 
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where we used the "eth" and "eth-bar" [12,1131, defined as: 



d d 

3 S = — [ "tt-t -\- i esc 6 — s cot 6 I = 8 + s c °t , (25) 

9 d 

3 S = — ( — — i esc — — h s cot 6 I = 3 — s c °t ^ , (26) 



89 dip 

in which the subindcx s indicates that they act on a quantity of spin weight s. These operators, inspired on the 
angular operators of quantum mechanics, were designed to act on the spin weighted spherical harmonics Y s ' m (9, ip). 
The spherical harmonics live on the sphere, where they form a functional basis, and are defined only for \s\ < I. When 
the 9 operator acts on the spin weighted spherical harmonic, it raises the spin weight of the harmonic: 

3 S Y s Lm = y/(l-s) (l + s + l)Y s+1 l - m , (27) 
and, when the 3 operator act on the spin weighted spherical harmonic, it lowers the spin weight of the harmonic: 



B s Y s l ' m = -V(l + s) (l-s + l)Y s ^ l ' m , (28) 

so the operator 3 is also called spin raising and 3 spin lowering. For a larger discussion on these harmonics and on 
how to use them to extract physical information carried by the gravitational wave see for instance [3 EH- For us, it 
is more suitable to work with the function $ = r ^4, which is expected to have a constant behavior in the regions far 
from the black hole due to the peeling theorem ( see @, Q), so that the final perturbation equation becomes: 

[□* +D flf | $ = 2A> 3 T 4 . (29) 

where there was only a change in the temporal radial operator, now taking the explicit form: 



d 2 . , , „ „ , d 2 . . . , d 2 „ _ _ d . „ ,„ . „ d . „M 

r 



□• = - (r 2 +2Mr)— + (r 2 - 2Mr) ^ + 4Mr ^ + 2 (2r + M) ^ + 2 (2r - M) ^ + 2- . 



This equation Eq. (|29|) , together with the definitions of the operators (|22H24[) . is the equation for the gravitational 
perturbation of a Schwarzschild black hole in Kerr-Schild coordinates due to an arbitrary test source T^. There are 
several interesting features worth noticing before proceeding further. First, that the equation (|30p is regular in all the 
domain of interest due to the choice of smooth penetrating horizon coordinates. Second, there is a complete splitting 
between the radial-temporal part and the angular one, with a principal part identical to the corresponding for the 
wave equation. This implies that the gravitational character of the equation is determined by the lower order terms 
(ie, first order derivatives and independent ones). The hypcrbolicity of the gravitational perturbation equation (|29|) is 
guarantied by the hypcrbolicity of the scalar wave, which was shown for instance in ref. [lj|. Finally, it is remarkable 
the simple form of the operators (|22H24[) . which have a different angular action on the corresponding projections of the 
perturbed source term. This suggests that, at least in some cases, there may be a preferred angular decomposition of 
such components. In the next section we will present a simple example that will show this feature. 



IV. DECOUPLING THE ANGULAR PART FOR DUST MATTER 



Let us consider the matter source to be described by a dust-like fluid, that is 

= pu^Uv , (31) 

where p is the rest mass density and the four velocity of the dust. Furthermore, we will consider that the fluid is 
infalling radially in the black hole, so this four velocity has only temporal and radial components, 

u» = (u°,u\0,0) , (32) 

which are functions of r and t only. The evolution of the fluid is described by the continuity equation for the current 
vector, J 1 * = pu^, and the conservation equation for the stress energy tensor: 

J M w - , (33) 

= . (34) 
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Althought the equations (|34|) imply in general the Euler equations, in the case of dust they simply reduce to the geodesic 
motion u u u^ v = 0. These geodesic equations, by using the symmetries of the spacetime and the normalization on 
the velocity = — 1, can be integrated once, so the four velocity can be expressed in terms of the constants of 

motion and the position of the particle (26[. For the radial motion that we are considering, we obtain the following 
expressions for the components of the four velocity: 



which determine the components of the velocity of each particle of the fluid at a given position and time. In these 
last equations E is a constant of motion along the trajectory, associated with the energy, and we will consider the 
minus sign for the u 1 component, as we will study a shell initially already infalling. Subsequently, the hydrodynamical 
problem reduces to the continuity equation ([55]) , which provides a evolution equation for the density, namely: 

9t (V^pu ) + driV^gpu 1 ) = • (36) 

where yf—g is the determinant of the four-dimensional metric. 

Let us now study more in detail the source terms appearing in the perturbation equation. Within the simple form 
of the four- velocity (f3"2"| . only the projections of the stress energy tensor along the light cone in the k-direction will 
not vanish, that is 

T km , = T m . m , = , T kk = (A;" u^f p = (u° + u 1 ) 2 p . (37) 



In this simplified case, the source terms for the perturbation equation (|29|) are just given by 



T 4 = T k k T k k = J u °+f) 2 Q_ ldop . (38) 

This simple expression indicates that the angular part can be decoupled in the source terms by decomposing the 
density in terms of the usual spherical harmonics (with zero weight), that is, 

p = Y,Ph™^r)Y l > m {6,<l>) ■ (39) 

lm 

The action of the bar eth operators (|26[) acting on an harmonic of spin zero lowers twice the spin weight to —2, 



B-^oYo 1 '" 1 = -VI (Z + l)3-iy-i'' m = VC-l) rO + 1) (Z + 2)y_ 2 '' m . (40) 
Collecting these results we get that the source term for this dust like case has the form 

T 4 = f kk T kk = - ( "°^ 1)2 X>.™(*,r) VG-l) 1(1 + 1) (l + 2)Y^ m . (41) 

Irn 

Notice that we are able to describe configurations of dust not necessarily spherically symmetric. The mass of the 
cloud will be given by the (I, m) = (0, 0) mode, which do not generate any gravitational perturbation, while the other 
modes will describe the over and under densities to this spheric base mode. 

The same decoupling of the angular part can be performed to the gravitational perturbation equation, by decom- 
posing the Weyl scalar <!> as a function of spherical harmonics with spin weight —2 0: 

$ = J2 R i,™(t,r)Y- 2 l > m (6,4>) , (42) 

lm 

such that the angular operator for the gravitational perturbation equation (|2ip is expressed in terms of the eth 
operators defined in Eqs. (|25|26j) in the following way: 

□ 0¥ , = cLig_ 2 . (43) 

Using the properties of the eth operators acting on the harmonics, given by eqs. ([27l [28]) , it can be seen that the 
Y_2 l,m are eigenfunctions of the angular operator: 

□ fl Y_ 2 l > m = §_! 3_ 2 Y_ 2 l ' m = — (£ — 1) (1 + 2) r- 2 '' m ■ (44) 
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Collecting these results we get that all the terms in the perturbation equation, including the sources, are multiples 
of Y_ 2 l ' m , so that we can obtain an equation for each mode (/, m) of the radial-temporal part of the 

- (r 2 + 2Mr) ^ + (r 2 -2Mr) +4Mr^ + 2 (2r + M) ^ + 2 (2r - M) + 

M \ / 'E-JE 2 - 1 + 2^ \ 2 , 

2— -(/-l) (Z + 2)J ^, m +4vrr 3 r _ 2M H V(2 -!)/(/ + !) (1+2)^ = 0. (45) 



An analogous expression can be found for the density modes from cq. (|36[) . where after substituting the decomposition 
(|39|) and the explicit expressions of the four- velocity (|35|) . one obtains 5 



E 2 -l 



3M 



d t pi, m + v r d r pt, m + 2 ^ — |^ r v r p;, m = , (46) 

where we have introduced the coordinate velocity, v r = ^? = ^ ■ It can be shown that all the terms are regular and 
well behaved in all the domain outside r = 0. Any initial arbitrary matter distribution can be expanded in terms 
of spherical harmonics with spin weight s = —2, although! neither the monopole nor the dipole modes will generate 
any gravitational reaction. This decomposition allow us to study, with a ID numerical code, any radially in-falling 
dust matter distribution and its gravitational reaction. Notice also that there is a degeneracy with respect to the m 
modes, since only the label I appears in the evolution equation (|45[) . This is due to the strong symmetry conditions; 
the background is spherically symmetric and we are restricting the fluid movement to be radially infalling. Indeed, 
the m modes give a description of the functions in the azimuthal angle tp, and due to the spherical symmetry of the 
background and of the motion of the fluid, we see that the action of the source on the black hole is independent of 
the azimuthal angle of the source. 

Thus, in this simple case of radial infall, each mode of matter awakes the corresponding mode in the gravitational 
radiation. In order to study the dependence of the gravitational wave on the initial distribution of the source which 
generates it, we will solve numerically the gravitational perturbation equation by recasting it as a system of first order 
evolution equations. In the next section we describe our numerical procedure to solve this system of equations and 
the results obtained from these evolutions. 



V. NUMERICAL EVOLUTION 



The Teukolsky equation (|45[) can be reduce to a first order system of evolution equations by following a standard 
procedure used for the scalar wave equation (l7j . which guaranties a well posed system of equations. We define the 
functions ipi t m an d II/. m in terms of the Ri, m functions as 

1pl,m = 9 r Rl,m , ^l.rn = [d t R^ m - (3 r 1pl m ) , (47) 

a 

so, using the lapse, shift and metric functions given in Eq. (|16|) for the Schwarzschild spacetime described in penetrating 
coordinates, we get for the definition of Ih jTO 

Ib. m = T + 2M d t Ri, m -2 — V>;, m , (48) 
r r 

and after a straightforward derivation for the evolution equation for n; m , derived from Eq. (|48l45p . we get the following 
first order system of equations: 

d t Ri, m = r M (rU 1 , m + 2Miln, m ) , (49) 

(1 \ 1 2 M 

r + 9M ( r n '> m + 2 M V'i.m) J = - - 2 M (r r Ul >m + 2 Md r 1pl tm ) + + 2 (Hf jTre - IpLrn) , (50) 

dtlli, m = l \ nr (2Md r Il Lm +rdr^ m )+ . , 2 an2 ((2r 2 + 5Afr + 4Af 2 ) n / . m + (r + 4M) (2r + 3M)Vi,m) 
r + 2M r(r + 2My 



E-JE 2 -1 + 2^. 



2^- a 1} / + 2) ) Ri, m + ^r I \_ 2M 
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This system will be evolved simultaneously with the continuity equation (|46[) for the density of the source by using 
the Method Of Lines (MoL) [l8| . Within this approach, the space derivatives are first discretized in order to obtain 
a semi-discrete system. The semi-discrete system is evolved in time by using (at least) a third order Rungc-Kutta. 
Due to the simplicity of our equations and the absence of shock, standard fourth order finite difference operators have 
been used for discretizing the space derivatives. In addition, a small Kreiss-Oliger dissipation has been added in order 
to damp the unphysical high-frequency modes. 

Our computational grid has some boundaries at finite positions. In general, data at those boundaries has to be 
provided in order to get a stable evolution. In our case, the boundaries are located at ri n = M and r out varies 
according to the case analyzed. The interior boundary ri n is placed inside the apparent horizon of the black hole. In 
this case, it is not necessary to impose a boundary condition, since all the characteristic fall into the black hole; all 
the information travels towards the same direction. The outer boundary r out is placed far away from the black hole, 
and we are imposing outgoing boundary conditions for the fluid density p and maximally dissipative BC for the <E>. 

We have consider as a initial data a shell of matter which can be written simply as the superposition of two spherical 
harmonic modes, namely 

p = p 0fi {r, t)Y °>° + p lim (r, t)Y l ' m , I > 2 . (52) 

The precise form of the / = mode po,o is not important for the radiation outcome, since only I > 2 modes of 
the density do contribute to the right-hand-side of the perturbation equation. In other words, the radiation is (at 
least) quadrupolar, so neither the monopole nor the dipole part produce gravitational radiation. Notice however that, 
although dynamically the (0, 0) mode is not important, it is the one that determines the total mass of the matter 
distribution. 

In this section, we start by presenting an study of a narrow Gaussian pulse of matter falling onto the black hole and 
generating a gravitational wave, and show that we obtain the expected behavior, validating in this way our equations 
and measuring the accuracy of our code. Next, we will consider three cases to study the gravitational response of 
the black hole to the infalling matter. In the first case we study a single pulse based on a gaussian distribution and 
will vary its width, describing the wave-form generated by the corresponding pulses and measuring its dependence 
on the width of the initial pulse. In the second case, we will describe three consecutive gaussian pulses and describe 
the gravitational radiation generated, in order to approach a model of matter falling onto the black hole in a periodic 
manner. Finally, we present a third case with two consecutive Gaussian pulses, and vary the separation between 
them, in order to study the waveforms of the gravitational wave as a function of the infall frequency. 



A. Testing the code 

We start by considering initially a single Gaussian pulse of matter: 

p l , m (r,t = 0)=A e- ( - r - r ^/° 2 , (53) 

where Aq,tq are the initial amplitude and the initial position of the center of the gaussian and a is its width. We 
will consider that the pulse is in rest at infinity, that is, E = 1. We do not consider the gravitational perturbation 
generated by the fact that the pulse arrive to the finite position rg, so at the initial time there is no gravitational 
radiation, that is R(t = 0, r) = ip(t = 0, r) = = 0, r) = 0. Due to this inconsistency, a spurious radiation is 
developed very quickly from the pulse and infalls to the black hole. The scattered part goes away to infinity, leaving 
the pulse of matter moving with a consistent set of data. The initial pulse of dust is located at vq = 70M, which is 
far enough in order to not interact strongly with the spacetime background until the initial spurious wave content left 
the domain. 

We assume first a density mode p2,o, with an amplitude Aq = 10~ 5 , and a width of a = 0.5. In fig. [T]is compared 
the profile of such mode at different times (i. e., every t = 25M) with the exact solution for the envelope of the 
density 

Pi, m (r)=A (^) § , (54) 

which is computed by integrating the evolution equation (|46[) for the mode pi, m - It can be seen that the fluid is 
evolving properly since the envelope fits during the evolution. In Fig. [2] it is displayed the radiation $ produced by 
this shell of fluid, measured by three different observers placed at r = 100, 150 and 200AZ . The signals has been 
shifted in order to overlap and stress the expected 1/r behaviour of the M/4 (ie, a constant behaviour of <E>). From the 
same signal we can adjust the ring-down form e~ Ul * cos {uj r t) in order to obtain the quasinormal frequencies {ujr, u>i}. 
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FIG. 1: The evolved density profile for a pulse of dust, plotted every t = 25M, and the exact envelope of the density. The 
gaussian has an initial amplitude of Aq = 10 -5 , an initial width of a = 0.5 and is located at ro = 70 M. 
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FIG. 2: Gravitational signal produced by the pulse from fig.[U measured by three different observers located at r = 100 M, 150 M 
and 200 M. The signals have been shifted in time in order to overlap them, as expected from the 1/r decay of the ^4. 



For this pulse, with an observer located at r b s = 100 M, we obtain the values M ui = 0.3734— 0.0891 i, which are 
very close to the analytical values Mui = 0.37367— 0.08896 i 0, [2(|. In fig. [3] it is displayed the same gravitational 
response $ in logarithmic scale in order to better appreciate the expected late time behavior, known as tail, which 
T7L 2l| . In fig. [4] is shown the radiation produced by the same initial parameters for the profile pi.m 
The corresponding QNM frequencies also are in good agreement with the theoretical ones 



goes as t 
for 1 = 2,3,4,5 



which are summarized in table I. 



B. One pulse with various widths 



Once we have tested the accuracy of our code by comparing with already known results, we proceed to analyze 
the gravitational response for different widths of the pulse. We studied the cases for a = 0.5 M,M, 1.5 M, 2 M, and 
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time 



FIG. 3: The $ from fig. [2] in logarithmic scale (continuous line) and a function proportional to t 7 (dashed line). Notice how 
the late time behavior of the gravitational response (ie, the tail) decays like t~ 7 . 



TABLE I: Quasinormal frequencies, notice that for the Schwarzchild case these values are independent of m 



I 


Computed Mui 


Analytical Mui 


2 


0.3734 - 0.0891 i 


0.37367 - 0.08896 i 


3 


0.5990 - 0.0932 i 


0.59944 - 0.09270 i 


4 


0.8055 - 0.0957 i 


0.80918 - 0.09416 i 



2.5 M, and present the corresponding logarithmic plot of the gravitational response in fig. El The outgoing energy 
flux per unit of time is given by: 



dE ,. r 2 
— — = hm 

dt 1 — >O0 16 7T 



^ d dt' 



dil = lim 

r^oo 16 7T 



Rl,m dt' 



(55) 



where we used the multipole expansion (|42|) for $ and the orthogonality of the spin weighted spherical harmonics. The 
total emitted energy can be computed by integrating in time the outgoing flux ([53)) . and the result for the different 
widths is displayed in fig. [5] . 




FIG. 4: Gravitational signal for different modes, namely I = 2, 3, 4. All of them have the same value for the initial amplitude, 
Aa = 10°, the same width, a = 0.5, and start at the same position, ro = 70 M. 




time time 

FIG. 5: The gravitational signal save for the different value of a. We are plotting, from left to right, for a = M, 1.5 M, 2 M, 2.5 M. 
Notice how the QNM frequencies become less noticeable, and the late time behavior differs from the theoretical one, as the 
matter pulse becomes wider. 

Our simulations show that the response of the black hole is the well know ring down only when the source of the 
perturbation acts in a very short period of time. As the source of the perturbation occurs during a longer interval of 
time, the signal loses the normal modes of the black hole and consequently the ring down behavior is gone. Also, with 
respect to the expected tail behavior, we see that the agreement with the theoretical exponential decay, gradually 
happens at later times as the matter signal is wider. These results have been also noticed by @, where they claim 
that "the QNM excitation is induced by the curvature profiles that have spatial wave lengths comparable to the width 
of the black hole potential" , which is a conclusion similar to our own. 

Indeed we see that in general, the usual gravitational signal has an initial burst which is then followed by a decaying 
oscillating wave with the ring down frequencies. Our simulations show that as the pulse becomes wider the signal only 
keeps the first initial bursts. This behavior might have important consequences in the actual detection of gravitational 
waves as long as the ring down is not present in all the gravitational emissions of a perturbed black hole, and should 
be confirmed by simulations performed with data resembling realistic astrophysical scenarios. 

C. 3 pulses 

Periodic and quasiperiodic variations are observed in different classes of astrophysical objects containing accretion 
discs. Quasiperiodic oscillations of thick accretion discs orbiting around a black hole have been addressed as sources 
of gravitational radiation [23| . In such cases, the tori oscillates induced by perturbations from the equilibrium 
configuration; as the inner ring of the tori approaches the black hole, part of its mass is moved through the cusp, 
resulting in a quasiperiodic accretion of matter onto the black hole. The dynamics of the tori is consequently printed 
in the gravitational wave signal. 

In order to emulate this scenario we set up a density profile composed by three gaussian pulses, mimicing a periodic 
perturbation like the one expected in the accreting tori [22| ■ Although we are not capturing the complete picture of 
the dynamics of the tori (like the angular momentum accreted) , this approach allow us to study the main features of 
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FIG. 6: Total energy as a function of the shell width, measured by an observer located at r = 70M. 

the gravitational signal produced by these periodic perturbations. 

We analyze the gravitational response due to three consecutive pulses, 

P2, (r, t = 0)=A (e-('-'-o) 2 /- 2 + e -(r-r -d)> 2 + e -(r-r Q -2df/a^ (56) 

with the closest one located at tq = 185M and the other separated a distance d = 15A/ behind each other. The 
initial amplitude is A$ = 1 x 10~ 5 and the width is a — 0.5M for each of them. We present in figure [7] (top panels) 
the gravitational signal measured by an observer located at r = 70M, where the periodic signal of the initial density 
profile is printed upon it. Notice however how the ring down behavior of the gravitational response is preserved if the 
pulses are far enough from each other. To make clear this point we have superposed the signal that is produced by 
a single pulse located at ro = 200M (ie, the central one) and notice that the general behavior coincides. This result 
differs from the one obtained by [22j in which they obtain that the gravitational response due to a periodic infall 
of matter generates a corresponding gravitational signal made of burst without the ring down behavior. This may 
happen when the infalling pulses arc so close to each other that there is no time for the appearance of the normal 
modes in the gravitational responses, as it is shown in fig [JJ (bottom panels) , where the separation between pulses is 
just d = 5M. This behaviour will be analyzed in more detail in the next subsection. 

D. 2 pulses 

In this section we consider the case of two consecutive pulses, varying the separation between them, in order to 
study in detail the transition between a ring-down response of the black hole and a forced oscillation due to the infall 
of quasipcriodic pulses. The initial data is a density profile of the form: 

P2, (r, t=0)=A (V(—o) 2 /- 2 + e ~(r-r -df/^ (5?) 

where d varies between [15 — 1]. The center of the first pulse is ro = 200 and we use a = 0.5M as usual. In figure 
[5] it is shown how the signal mirrors the behavior of the shell that is falling into the black hole; as the separation 
between the pulses decreases, the time elapsed between the signals measurement also decreases, up to the point where 
the observer is unable to distingish between two consecutive signals. This happens in the second and third panel of 
the plot, where each signal has lost their individuality superposing each other, to form a signal that hides the ring 
down behavior of the individual. Finally in the fourth panel the pulses are initially so close 1M that the gravitational 
response of the hole is the same as produced by a single perturbation with the double of the width. Based on these 
results we can conclude that the initial separation of the matter pulses actually affects the gravitational response of 
the black hole. In particular, the signal corresponding to the pulse arriving first when the second one is close (5M) 
has no time to generate a ring down signal. The last pulse however does present the ring down behavior as expected. 
In particular when the initial separation is 5 M the signal resembles a succession of bursts similar to the ones obtained 
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FIG. 7: Gravitational response due to the infall of three consecutive pulses. In the top panel we are showing the $ for three 
pulses with a separation d = 15M, superposed with the response of one single pulse located at r — 200Af in the left panel. In 
the bottom panels the separation between pulses has been reduced to d = 5M. 

by [22j . Finally, we notice that when the pulses are even closer to each other the gravitational response becomes the 
one due to a single wider pulse. Regarding the separation of consecutive pulses, there is a similar behavior as in the 
case of the width of the pulse. The gravitational signal is affected when the initial separation between the pulses is 
comparable to the radius of the black hole. 



VI. DISCUSSION 



We have rederived the evolution equation for gravitational perturbations including sources, showing how to deal 
in a consistent manner with different choices of signature. As a first application, we have studied the spherically 
symmetric static spacctimcs, in particular the Schwarzschild black hole. After some comparison with previous similar 
works, we have constructed the tetrad and choose the penetrating coordinates to compute the perturbation equation 
with matter source terms. 

By means of the expansion of the rest mass density in spherical harmonics, we were able to describe a generic 
distribution of dust radially falling onto the black hole as a one dimensional problem. We have implemented the 
perturbation and the matter equations numerically to study the gravitational response in several interesting cases 
of the infalling dust. We have first discussed the importance of the width in the infalling shells in the generated 
gravitational signal, showing how the ringdown behavior is lost in favor of a single burst response, when the pulse 
becomes wider. Afterwards we have studied the case of infall of consecutive pulses. Our claim is that this is a way 
to study the gravitational perturbation of a black hole, when the matter is released in a periodic fashion, as in an 
oscillating tori. Indeed, we are thinking on a tori which initially fills a region just down to the Roche's lobe and 
due to some perturbation, starts to oscillate and in each period, some of its matter gets further than the Roche's 
lobe, consequently falling onto the black hole and producing a periodic source of gravitational perturbation. We 
have modeled the periodic emission of matter by the tori as consecutive pulses of radially infalling dust. Moreover, 
once the matter passed beyond the Roche's lobe, its motion could be consider is practically radial, so our approach 
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FIG. 8: It is shown the gravitational response of two consecutive pulses varying the initial separation d between them. These 
separations are, from left to right an up down, d — 15M, d — 5M, d — AM and d — IM. 



allow us to focus in this stage of the motion. There were claims that in these cases of periodic infall of matter, the 
gravitational response was modified losing the ring down behavior in favor of a response dominated by the frequency 
of the infalling matter as in a forced motion, we started with three consecutive pulses and obtained a gravitational 
signal with three burst each one with the spected ring down behavior, with no further effect of the periodicity of 
the infalling matter. Finally, we proceded to study the gravitational response of two consecutive infalling pulses 
varying the initial separation and obtained that indeed there is an interval of the separation of the pulses in which 
the gravitational response is certainly affected, clearly when there is no time for the ring down to happend due to the 
infall of the next pulse, however if the initial separation gets even shorter the gravitational response is that of a single 
wider pulse. In this way we have showed that the gravitational signal is not only determined by the parameters of 
the space time it moves but it has also imprints of the manner in which it has been generated 

The present work gives a solid base and sets us in a track that should be follow on other cases, like the one with 
matter or a small compact body falling onto a rotating black hole. The presence of certain types of pressure can 
give rise to oscillating motions of spherical shells, as discussed in [HI, HH, which could generate a continue source of 
perturbation to the black hole, which in turn might generate a continued emission of gravitational waves, carrying 
information of the background as well as of the matter which produced them. 
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